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Abstract 

A mystery surrounds the stability properties of the splay-phase periodic solutions 
to a series array of ./V Josephson junction oscillators. Contrary to what one would 
expect from dynamical systems theory, the splay state appears to be neutrally stable 
for a wide range of system parameters. It has been explained why the splay state 
must be neutrally stable when the Stewart-McCumber parameter (3 (a measure of the 
junction internal capacitance) is zero. In this paper we complete the explanation of 
the apparent neutral stability; we show that the splay state is typically hyperbolic - 
either asymptotically stable or unstable — when (3 > 0. We conclude that there is only 
a single unit Floquet multiplier, based on accurate and systematic computations of the 
Floquet multipliers for (3 ranging from to 10. However, N — 2 multipliers are extremely 
close to 1 for (3 larger than about 1. In addition, two more Floquet multipliers approach 
1 as (3 becomes large. We visualize the global dynamics responsible for these nearly 
degenerate multipliers, and then estimate them accurately by a multiple time-scale 
analysis. For N = 4 junctions the analysis also predicts that the system converges 
toward either the in-phase state, the splay state, or two clusters of two oscillators, 
depending on the parameters. 

Abbreviated Title: Stability of periodic solutions in Josephson arrays. 
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1 Introduction 



Arrays of Josephson junctions have been studied intensely by mathematicians and physicists 
because of their interesting properties as dynamical systems. A single Josephson junction 
oscillator is analogous to the familiar physical pendulum, and the coupling of an arbitrary 
number of oscillators allows us to study the dynamics of the extended systems and the 
connection between ordinary differential equations and partial differential equations. On 
the other hand, there is also a large practical interest in Josephson junctions, which stems 
from their extremely high frequencies and sensitivity to the magnetic field | |J JA88| , |ASC94| . 
The typical oscillation amplitude of a single junction is so small that many junctions must 



be coupled together in order to get a macroscopic voltage oscillation |Jaifcal84j , [Hadfcal88 
BenfcBur91| 



One-dimensional arrays of Josephson junctions may be classified as either series or par- 
allel. Parallel arrays are well approximated by a nearest-neighbor diffusive coupling. This 
leads to discrete sine-Gordon models, which support "solitons" propagating across the junc- 
tions (see [|Watfcal95|1 and references therein). In series arrays, the junctions are coupled in 
an all-to-all fashion. For this "global" coupling, the spatial coordinate is irrelevant, unlike 
in parallel arrays, and particularly good analytical progress has been made. 

In this paper, we consider series arrays of N identical Josephson junctions shunted by a 
series LRC load (Fig. [I]). After standard non-dimensionalization, the governing equations 
reduce to [ |Hadfcal8l3| : 

f3(pj + (f>j + sin fa + Q = I b , j = 1, . . . , N, 

LQ + RQ + ^ = ^J24> k . (1) 
° iv k=l 

Here, <pj is the phase difference across the j-th junction, lb is the DC bias current, and NL, 
NR, C/N are the load inductance, resistance, and capacitance, respectively. (Thus L is 
the inductance per junction, etc. This scaling, introduced in |Hadfeal8"B| , is convenient for 
comparing different numbers of junctions.) As mentioned in the previous paragraph, the 
junctions are globally coupled even though they are arranged in series. Mathematically, the 
equations have permutation symmetry with respect to the indices j. The only nonlinearity 
is the sin term. In particular, the second equation is linear. We emphasize that the 
bias current lb is constant in time. Josephson junction arrays are often driven by such a DC 
forcing since it can still generate useful AC responses. In addition the equations become more 
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Figure 1: Circuit diagram for a series array of N identical Josephson junctions subject to 
a constant bias current h, and shunted by a series LRC load. 

difficult to analyze with an AC bias current. (The US voltage standard uses a series array 
of Josephson junctions with an AC bias current: The average voltage across the junctions is 
proportional to the frequency of the bias current due to the Josephson relation ||Kau&;al8"7|| .) 
The dimensionless Stewart-McCumber parameter /3 is a measure of the junction internal 



capacitance. Junctions can be made with a wide range of j3, ranging from 10 to 10 ||JJA88 



ASC94] , vdZ&;al94| . Much progress has been made in analyzing the system ([!]) in the limit 
(3 = 0, where the phase space of each oscillator is a circle ||Jaifeal84] , |Tsafeal9T| , |Swifeal92 



[rsafcSch"9^ , |StrfeMir93| , |WatfcStr94|| . In this paper we study the more difficult case when 
j3 > 0, where the shape of the oscillation in the (<pj, <fij) plane is affected by the coupling. 

A mechanical analog of a single Josephson junction (Eq.(^) without Q) is a damped 
pendulum of mass (3, with applied torque h- The coupling in Eq.([l]) arises in the mechanical 
analog if the pendula slide with friction around a common axle, and the axle is itself a 
torsional harmonic oscillator (Q is the angle of the axle, L is proportional to the moment of 
inertia, etc.). 

In this paper we consider only periodic oscillations that are overturning, i.e. (f)j(t + T) = 
<pj(t) + 2tt for all t and all j where T > is the period. The most important periodic 
solutions of system (|l|) are the in-phase and the splay-phase states. The junctions oscillate 
coherently in the in-phase state: <f>j(t) = 4>\{t) for all j. The splay state is also periodic, but 
the individual waveforms are staggered equally in time: <f>j(t) = (f>\{t+ (j — 1)T/N) for all j. 
(This is one splay state; the indices j can be exchanged arbitrarily due to the permutation 
symmetry of Eq.(|T]).) The in-phase state, due to its obvious coherence, has been considered 
promising as a high-frequency oscillator | |Lik86| , |Hadfeal8"%| , |MeyfeKre9U| ] . On the other hand 
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the splay state emits little rf power output, but other applications are sought, e.g. as a voltage 
amplifier |[lerfcBea95| . Since the existence of a single splay state automatically implies the 
existence of (iV— 1)! splay states by the permutation symmetry, the possibility of using these 
states as multiple- memory device has also been considered [ 0ts91 , (5chfeTsa9~i . 

The linear stability of periodic solutions is quantified in terms of Floquet multipliers. A 
periodic solution of an ordinary differential equation (ODE) in n dimensions has n Floquet 
multipliers. One of these multipliers is always 1, corresponding to the trivial perturbation 
along the solution, and the other n — 1 multipliers are the eigenvalues of the fixed point 
of the Poincare map. (The periodic solution of the ODE is a fixed point of the Poincare 
map.) A Floquet multiplier inside the unit circle in the complex plane indicates a decaying 
perturbation, while those outside the unit circle are "unstable." 

The linear stability of the in-phase state was numerically studied ||HadfcBea87| , |Hadfcal88| 
for various combination of the loads, the bias current If,, and (3. The results were explained 
analytically in the limit (3 = ||Jaifcal84| , |Lik86|| , then for a finite (3 > ||MeyfcKre90| , 
|(JhefcSch95|| using perturbation methods. The multipliers of the splay state were also studied 
intensively, as we review in the following. Roughly stated, the linear stability properties of 
the two solutions appear complementary; for example, when the in-phase state is stable, 
the splay state is usually unstable (although bistable situations or chaotic attractors may 
sometimes arise). 

However, many numerical studies made a puzzling observation: often, the splay state 
is not asymptotically stable when the in-phase state is unstable; rather, it is only neutrally 
stable (i.e. Liapunov stable but not asymptotically stable). For an .R-load with (3 = 0, Tsang 
et al. [[Tsafcal91|l found evidence that there are N Floquet multiplers on the unit circle (one 
complex conjugate pair and N — 2 unit multipliers). Tsang and Schwartz found that there 
are N — 2 unit multipliers for the (3 = junctions coupled with an LC-load ||Tsa&;Sch92[| . 
Nichols and Wiesenfeld ||NicfcWie92|| did an excellent study of five different cases: two with 
(3 = and three with f3 > 0. They always found, to within their numerical accuracy, N — 2 
unit multipliers when (3 = 0. For (3 > 0, they found that the splay state seems to be 
neutrally stable in iV — 2 directions for an LC or an LRC-lo&d, but asymptotically stable for 
a C-load. This difference is surprising because the C-load is a limiting case of the LC-load. 
Nichols and Wiesenfeld were very careful when stating their conclusions. They warned that 
numerical evidence cannot prove neutral stability because one cannot distinguish between a 
unit multiplier and one that is extremely close to 1. 
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All previous analytical studies of this peculiar phenomenon have assumed j3 = for 
simplicity. The neutral stability when (3 = can be understood in the weak coupling limit 
||Swifeal9"2] , |WiefeSwi95|| or in the limit iV — > oo ||StrfeMir93|| . A rigorous explanation was 
obtained in | WatfeStr94 |, where it was shown that the phase space of the system (|l]) has 
a simple foliated structure when (3 = 0. This structure gives rise to N — 3 constraints 
(constants of motion) to the system, and hence imposes the splay solution to have N — 2 
unit multipliers. This global result not only explained the peculiar linear stability of the splay 
state, but also made it possible to analyze the asymptotic behavior from an arbitrary initial 
condition. (Because initial conditions of the junctions cannot be specified, it is important 
that a desired periodic solution is globally stable, not just linearly stable, for any applications 
|Bra&Wie94j| .) In the weak coupling limit |WatfeStr94j| it was shown that the possible 
asymptotic states are either the in-phase state or one of the "incoherent states" which include 
the splay state special case. 

The present paper studies, numerically and then analytically, what happens to the phase 
space structure when (3 is positive, and connects the global picture to the local stability 
properties of the in-phase and the splay periodic solutions. Since Nichols and Wiesenfeld 
found apparent neutral stability in two of the three cases when (3 > 0, an early conjecture 
was that the simple structure of the (3 = equations might persist for (3 > 0, and one might 
hope to prove that the splay state is still neutrally stable. 

But as we show numerically in Section 2, the splay solution is not neutrally stable when 
(3 > with an LC-load. Unlike Nichols and Wiesenfeld |[NicfcWie92|| , we choose L = 1 
and C = 1/4 since this is "closer" to a C-load (1/C is large compared to L). We use 
N = 4 junctions since this is the smallest N that shows the degenerate structure of the 
(3 = equations. We visualize the global dynamics by using the time delay technique, and 
find that the splay solution is asymptotically stable for some set of (3 and This leads 
us to suspect that the neutral stability observed when (3 > in other regions of the (L, C) 
parameter space is also only approximate. Furthermore, an LRC-loa.d is qualitatively similar 
to an LC-load; the system is dissipative even with an LC-load due to the internal resistance 
of the junctions. 

The visualization of the flow in Sec. 2 also reveals that there are several different time 
scales in the dynamics when the splay state is stable. A trajectory initially converges to a 
manifold of "incoherent states", then slowly drift toward the splay state. This process can 
be quantified by monitoring the "moments" R m of the phases (defined in Sec. 2.2). The 
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first moment R\ converges to rapidly, then the second moment R2 slowly drift to 0. By 
choosing another set of (3 and /&, we can change the direction of the drift, while still forcing 
the incoherent manifold to be stable. R\ still converges to 0, but then R2 drifts toward 1. 
The final state is a periodic state with 2 antiphase pairs of oscillators. For yet another (3 
and R\ converges to 1, in which case the in-phase state is attracting. 

For (3 = 0, the drift of R2 is prohibited due to the rigid foliation structure; the drift 
is an important effect of positive (3. This change in the global dynamics is reflected in the 
Floquet multipliers around the splay state. Generically, there should be a single unit Floquet 
multiplier, but there are N — 2 unit multipliers when (3 = 0. For (3 > 0, N — 3 of them 
are perturbed away from the unity, and the genericity is recovered. With N = 4, only one 
multiplier is perturbed, which we call the "critical multiplier." 

To test this qualitative picture more quantitatively, we present in section 3 systematic 
computations of the Floquet multipliers as a function of (3. We choose N = 4 and the 
same load values as used by Nichols and Wiesenfeld |[NicfcWie92|| for both the LC and C- 
load. We continue the branch of splay solutions between (3 = and (3 = 10 using AUTO86 
poefcKer86|1 . This allows us to track the single critical multiplier into the regime where the 



splay solution is unstable. We find that this multiplier approaches 1 as (3 — > 0, as expected. 
However this multipler is also asymptotic to 1 as (3 — » 00. In fact, the multiplier approaches 
1 so quickly that Nichols and Wiesenfeld could not distinguish it from 1 when (3 = 1.1 in the 
LC-load case. (Our calculations are consistent with theirs.) At smaller (3, where the critical 
multiplier is significantly smaller than 1, the splay state is unstable in a different direction 
and thus Nichols and Wiesenfeld did not find the state with their forward integration. On 
the other hand, in the C-load case it turns out that the splay state is stable for (3 small, 
and they observed an asymptotically stable splay state at (3 = 0.1. For both loads, we also 
find that iV Floquet multipliers approach 1 as (3 — > 00, indicating that the junctions become 
uncoupled in this limit. 

There are two limits in which the oscillators become uncoupled regardless of the load: 
(3 — > 00, or lb — > 00. (Furthermore, the coupling is weak when the load impedance becomes 
large, for example R — > 00.) We follow [|ChefeSch95|1 , and consider the limit h — > 00. In 
section 4 we use the perturbation parameter e = 1/4 <C 1 to analytically verify the global 
picture visualized in Sec. 2, and estimate the Floquet multipliers computed in Sec. 3. Based 
on our observation in Sec. 2, we introduce three different time scales, and derive the order 
e A averaged equations for the phase drift of the oscillators. 
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In Sec. 4.3 the order e 2 averaged equation is obtained, which turns out to have the 
same form as the averaged equation when (3 = ||WiefeSwi95||: it predicts that the splay 



solution has N — 2 unit multipliers. While for (3 = this approximate result can be made 



rigorous [ WatfeStr94 |, we expect higher order corrections to all but one unit multiplier. The 
0(e 2 ) averaged equation controls the flow between the in-phase state and the incoherent 
manifold. The first moment Ri is shown to converge toward or 1, depending on parameters. 
The stability of the splay solution to perturbations of R\ corresponds to a single complex 
conjugate pair of Floquet multipliers. We estimate this multiplier pair from the second order 
averaged equation and find excellent agreement with our numerical results around the splay 
state. 

In order to resolve the drift along the incoherent manifold, we expand to higher order 
and derive the order e 4 corrections to averaged equation in Sec. 4.4. The second moment 
i?2 n ow converges to or 1. We compute the critical multiplier for N = 4, and again find 
excellent agreement with our numerical data. 

A summary and a discussion are presented in Sec. 5. Appendices A and B show details 
of the analysis in Sec. 4. 

2 Numerical integration of the equations 

In this section we numerically integrate the equations ([!]) using a standard ODE integrator. 
This has the disadvantage that we cannot find unstable states. Nonetheless we will convince 
the reader that the degeneracy observed when (3 = does not persist when (3 > 0, although 
the degeneracy re-appears approximately when (3 is large. 

Since Nichols and Wiesenfeld [|NicfcWie92|1 observed an asymptotically stable splay state 



in the C-load case, we are interested primarily in the LC and LRC loads. These are not 
qualitatively different in the Josephson junction system because of the internal resistance of 
the junctions. Thus, we only consider the LC load. 

We shall find that in some parameter regions of the LC-load system, there are asymp- 
totically stable splay solutions, which persist when a small nonzero resistance is added to 
the load. We will also show other types of periodic solutions that can be attracting. 
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2.1 Time-delayed coordinates 

Both the in-phase and the splay states are phase-locked states, and the phase relations among 
junctions are defined in terms of time lags rather than mere differences in angles <pj at each 
instant t since the waveform 4>j(t) is nonuniform. A process to convert a solution <pj(t) into 
the phase information is therefore necessary. 

One useful algorithm ||AshfcSwi93|| is based on time delays between the "firings" of os- 



cillators. We define each junction to "fire" when <pj = 0* (mod 2tt) for some fixed <fi* . (We 
choose <p* = 0.) Let tj(p) be the time of the p-th firing of the j-th junction. We monitor 
firings of all the junctions, and calculate the phases every time junction N fires. Assume 
that p and q satisfy tj(q) < tjv(p) < tj(q + 1). Then we define the phase (normalized time 
lag) of the j-th junction to be: 



(2) 



where Tj is an approximate period of the junction. The phase 0j(p) is approximately the 
angle that oscillator j covers between its most recent firing and the pth firing of oscillator N. 
Among the candidates for Tj, we choose Tj = tj(q + 1) — tj(q) in the following. This has an 
advantage that 9j always falls in the interval < 9j < 2n. The "time-delayed coordinate" 
9j is an angle coordinate and we identify 9j = with 9j = 2n. 

2.2 Periodic attractors 

We integrate the system ([!]) with iV = 4 by an adaptive time-step Runge-Kutta integrator, 
and compute the time-delayed coordinates. We have tested both LC and C loads, but will 
focus on one representative case when L = 1, R = 0, C = 1/4. An array with this load 
exhibits a variety of periodic solutions depending on the parameters (1^0). The initial 
conditions are chosen to illustrate the transient dynamics. 

Four cases are shown in Fig. |^. For each case, two graphs are shown; an upper graph 
shows the time-delayed coordinates 9j. Since 9^ = by definition, only three curves can 
be seen. Each lower graph shows the amplitudes R\ and R2 of the the first and the second 
(complex) moments z m of 9j's: 

1 N ^ 

Zm ^ t 6 k 1 Rjn \Zm\ ■ (3) 

k=l 
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Figure 2: Time-delayed coordinates 6j (upper graphs) and the amplitude of the first and the 
second moments R\ and R2 (lower graphs) as functions of t, for N = 4 Josephson junctions. 
Cases (a-c) uses L = 1, R = 0, C — 1/4, and If, = 1.5. Only the McCumber parameter (3 is 
changed: (a) 0.7, (b) 0.15, (c) 0.3. Case (d) uses the same load but with I b = 2.5 and (3 = 1. 
A trajectory converges to (a) the in-phase state, (b) the splay state, (c) the 2+2 incoherent 
state, and to (d) a (generic) incoherent state. 

The amplitude of the first moment, R±, was originally introduced by Kuramoto as an "order 



parameter" to measure the phase coherence of a population of phase oscillators [ |Kur84 
Indeed, R± = 1 if and only if the oscillators are in-phase. When R\ = 0, we say that the 
oscillators are "incoherent" (this is our definition of an incoherent state). The splay state is 
characterized by R m = 1 if m is a multiple of N, and R m = otherwise. The horizontal t-axis 
for the figures denotes the (discrete) firing time t^ip)- The first several firings are ignored 
since the junctions rapidly change their phases. During this transient, the time-delayed 
coordinates are not meaningful. 

For Figs. 0(a) to (c), the bias current is fixed at If, = 1.5, and only the McCumber 
parameter (3 is changed. Fig. 0(a), which shows a stable in-phase attractor, uses (3 = 0.7. 
Since 9j = is identified with dj = 2tt, we see that four junctions have become synchronized 
after t ~ 100. This is also clear from the moments: both R± and R2 converge to 1. 
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In Fig. 0(b), (3 = 0.15 is used. The phase differences converge to 2n/N, thus the attractor 
is the splay state. The way the trajectory approaches the attractor is somewhat interesting. 
Note that R\ vanishes at t ~ 200, after which the state is incoherent. If the angles 9j are 
plotted on the unit circle, they are scattered so that their center of mass lies at the origin. 
For N = 4, an incoherent state has two "antiphase" pairs: The oscillators can be ordered so 
that 

(01, 02, 03, 04) = (*!, Ml +M) (4) 

with < 0i < Ti. After Ri vanishes, only the angle 0i (and hence 03) in equation (|) 
changes, and the incoherent state is characterized by i?2 = | cos 6\ \ . Note the slow drift of 
i?2^) toward 0, indicating the splay state with Q\ = tt/2. 

Fig. 0(c) uses an intermediate value (3 = 0.3. At t ~ 200, Ri vanishes and the phases are 
again incoherent. A slow drift of R2 follows, as in (b) occurs, but in the opposite direction. 
At t ~ 300, i?2 converges to unity, indicating that Q\ approaches or tx. The final state 
is two clusters located 180 degrees apart. We call this state the "2+2 incoherent" periodic 
solution. 

With I b = 1.5 as in Figs. 0(a-c), we did a parameter sweep of A/5 = 0.01, and found that 
the splay state is stable for (3 < 0.21, the 2+2 incoherent state is stable for 0.19 < (3 < 0.61, 
and the in-phase state is stable for (3 > 0.40. This bistability between the splay and 2+2 
incoherent states implies that there exists an unstable generic incoherent state when 0.19 < 
(3 < 0.21. By "generic," we mean that < Q\ < tt/2 in equation (|j). 

Fig. 0(d) uses the same load, but both I\> = 2.5 and (3 = 1 are larger than the previous 
values. We see that R\ vanishes at t w 300 again while R2 stays constant for t > 250. The 
constant value is neither nor 1 so that the final state appears to be a generic incoherent 
state with < 0i < n/2. Different initial conditions converge to incoherent states with 
different R2 (and thus different X .) It appears that there is an incoherent periodic state for 
each value or R 2 (or 0x), just as there is in the equations with (3 = 0. However, we believe 
that the states are not exactly periodic but are drifting very slowly towards the splay state. 
In Sections 3 and 4 we will give evidence for this slow drift. 

When lb or (3 is larger, for example h = 2.5, (3 = 4 or If, = 4, (3 = 1, there appears 
to be a periodic state with any values of the angles 0j. During the transient dynamics the 
delay coordinates are not meaningful, but after this the angles 8j stay constant at values 
determined by the initial conditions. In other words, the oscillators are effectively decoupled, 



10 



and they maintain any phase relationship. Mathematically speaking, there appears to be 
an attracting 4-torus foliated by periodic orbits. As with the case described in figure 0(d), 
however, we expect that the dynamics are actually "generic" , but the drift of both R\ and 
i?2 occurs on a time scale too long to observe. 

All of the oscillators have the same waveform in the in-phase, splay, and 2+2 incoherent 
solutions. The existence of such solutions has been proved in [|ArofcHua95|| for a wide range 



of loads, including ours. In a generic incoherent state implied above, however, one antiphase 
pair has a different waveform from the other pair (see section 2.4) and there is no existence 
proof for this type of periodic solutions. 

2.3 Three-dimensional projection 

The dynamics become more apparent by using a linear transformation of the time-delayed 
coordinates ||AshfeSwi92| , |TsafeSch92| , |WatfeStr94|| . Since we are only interested in the rela- 



tive phase differences, we change the coordinates by, for example, m ,i = Q\ + 9 2 ± 3 ± #4 and 
^2,3 = ±#1 + #2 + #3 — 84, then only show u\, 112, and M3. Figs. |||(a-d) show the dynamics of 
Figs. §(a-d), respectively, in these coordinates. 

The trajectories run inside a tetrahedron in this representation. The faces of the tetra- 
hedron correspond to when two 9j's are identical (a "2+1+1" state). Each edge of the 
tetrahedron corresponds to one of two types of oscillation: "2+2" (shown as dotted lines) 
or "3+1" (as solid lines). Since no change in ordering of junctions occurs in Fig. ^| after a 
short transient, the phases maintain the relationship, say, 64 = < 63 < 62 < Q\ < 2n. 
Thus, a trajectory does not escape from the tetrahedron. (Such an escape is prohibited if 
there is an attracting invariant 4-torus. However, with strong coupling the oscillators can 
shuffle their order even after they have converged to the attractor.) The four vertices of 
the tetrahedron correspond to the in-phase state. (There are four ways to approach this 
state while maintaining the ordering.) The incoherent bar inside the tetrahedron is the set 
of incoherent states for which Ri vanishes. (The incoherent bar, indicated by a dashed line, 
generalizes to an (N — 3) -dimensional incoherent manifold when there are N oscillators.) 
The two end points of the bar are 2+2 incoherent states, while the dot at the midpoint is 
the splay state. The trajectories shown are actually discrete sets of points at {^(p)}, but 
the points are connected in the figures for clarity. 

In Fig. |3|(a) the trajectory starts at a point inside the tetrahedron, but converges to the 
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Figure 3: Projection of trajectories onto a three-dimensional space denned by a linear 
transformation of the time-delay coordinates: tt ,i = Q\ + 2 ± 6*3 ± 6* 4 and w 2 ,3 = ±#i T #2 + 
#3 — #4- Only the phase differences U\, u 2 , u 3 are shown. The cases (a-d) correspond to the 
cases in Fig. |2], respectively. The vertices of the tetrahedron correspond to the in-phase state 
while the bar inside the region is the set of incoherent states. The two endpoints of the bar 
are both 2+2 incoherent states, while the dot in the middle of the bar is the splay state. 



in-phase state at a vertex. In Fig. 0(b), on the other hand, the initial condition is near the 
in-phase state but the trajectory converges to the splay state. This process consists of two 
stages; the first fairly rapid fall onto the incoherent bar is followed by a slow drift along the 
bar. 

In Fig. 0(c) the initial condition is also near the in-phase state, and the trajectory falls 
onto the incoherent bar. However, this time, it converges to the 2+2 incoherent state at the 
end of the bar. In Fig. 0(d) the trajectory does not appear to move along the incoherent bar 
once it hits the bar. The calculation is halted at t ~ 600. However, we expect that there is 
a very slow drift along the bar as described earlier. 

The drift along the incoherent bar implies that the degenerate phase space structure of 
the P = limit is broken for (5 > 0. In this limit, the tetrahedron in Fig. is foliated by a 
family of invariant subspaces |[WatfcStr94|| . If such degeneracy were to persist when (3 > 0, 
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then the structure would prohibit a trajectory from moving from one subspace to another, 
thus the drift would not occur. As a result of this drift, the splay state, the 2+2 incoherent 
state, or more generic incoherent states can become asymptotically stable. 

Is the incoherent manifold dynamically invariant? For N = 3 the incoherent manifold 
is precisely the splay state, so the answer is yes. Appendix B gives a formal proof that the 
incoherent manifold is dynamically invariant for N = 4 and weak coupling. The proof is not 
rigorous because it relies on the phase-shift symmetry induced by averaging. In Appendix 
B we also show that the incoherent manifold, defined by R\ = 0, is not invariant for N > A. 
However, when the coupling is weak there is an invariant manifold, close to the set Ri = 0, 
on which there is a slow drift. Thus, the concept of an incoherent manifold is still useful. 

2.4 Waveforms 

The time-delayed coordinates 9j are useful to investigate the dynamics, but do not show di- 
rect information on <f)j(t). Therefore, we investigate the waveforms and symmetry properties 
of the periodic solutions in this section. 

Figures |](a-d) show the waveforms for (a) the in-phase state, (b) the splay state, (c) the 
2+2 incoherent state, and (d) the generic incoherent state, respectively. (The case (d) is 
close enough to being periodic.) Each upper graph displays sin </>.,■ (the supercurrent through 
the j-th junction) as a function of t, whereas each lower graph shows the same solution in 
terms of the phase portrait (0j,0j) for all j. The dashed curve in the lower graph denotes 
the trajectory of an uncoupled (single) junction with the same parameters (Ib,f3)- 

The trajectories for the in-phase state and the 2+2 incoherent state are significantly 
altered from the trajectory without coupling. In contrast, the trajectories for the splay state 
and the generic incoherent state are almost unchanged in shape. This is partly due to the 
parameter values; note that figures (b) and (d) have the largest j3 values. Furthermore, when 
the state is incoherent the coupling terms in equation ([I]), J2(pk, tend to cancel each other. 

In all four cases, the four trajectories in the phase portraits follow the same curve. This 
must be the case for all but the generic incoherent state in figure (d). Let us explain: We have 
already given the general form of the in-phase, splay, and 2+2 incoherent solutions. Each of 
these states has a single waveform, with different phase shifts between the oscillators. For any 
incoherent state with N = 4, the oscillators can be numbered so that 4>3(t) = 0i(t+T/2), and 
<pA_{t) = 4>2{t + T/2), where T is the period. In figure (d), we observe that (p2{t) = 4>i{t + c) 
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Figure 4: The supercurrent sin0j(t) (upper graphs) and the phase portrait <pj versus <pj 
(lower graphs) of the final state achieved in Fig. |2], respectively. The dashed curve in the 
phase portrait shows the trajectory of an uncoupled junction with the same lb and (3 values. 

for some phase shift c, but in general these two waveforms can be different. We have in 
fact observed two distinct waveforms for a generic incoherent state in a different parameter 
regime. 

Finally, we show a way to distinguish between the different types of oscillations which 
does not require the time-delay technique. In Fig. || we plot 2 — 04 vs. 4>i ~ 03- Since (f)j is 
proportional to the voltage across the Josephson junction, we call this the Lissajous figure 
created by the voltage signals. 

All <pj(t) are synchronized for the in-phase solution, thus the figure for case (a) appears as 
the point at the origin. The pattern of the other states depends on the order of the oscillators: 
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Figure 5: Lissajous figures obtained by plotting the difference of the voltages across junc- 
tions 2 and 4 versus that of junctions 1 and 3. Four figures are superimposed, corresponding 
to the periodic states in Fig. |2|. The in-phase state (a) appears as a point at the origin. 

we choose initial conditions or permute the oscillators so that <f>\ < 02 < 03 < 04 < 0i + 27r. 
The splay state (b) has a 4-fold symmetry while a generic incoherent state (d) has only 2-fold 
rotational symmetry about the origin. The 2+2 incoherent state appears as a diagonal line 
with a slope +1 (The slope can also be —1, depending on initial conditions.) The symmetries 
of the Lissajous figures are easily explained from the waveforms. For example, phase shifting 
the splay state by a quarter period rotates the Lissajous figure by 90 degrees. 

The Lissajous figures can be extended to three dimensions using the u z" coordinate 
0i — 02 + 03 — 04- The resulting figures are analogous to Fig. ||, but they can be made 
without computing the time-delay coordinates. The z coordinate of the Lissajous figure is 
needed to distinguish between an incoherent oscillation (on the incoherent bar in Fig. and 
a general oscillation (in the interior of the tetrahedron). Both oscillations look like ellipses 
in the 02 — 04 vs. 0i — 03 plane. However, incoherent oscillations have the extra symmetry 
that 0i — 02 + 03 — 04 is unchanged after half a period. 

3 Floquet multiplier computation 

In this section we report on calculations of the Floquet multipliers of the splay-phase state. 
The Floquet multipliers provide local information which is complementary to the images 
of the global phase space shown in the previous section. The degenerate dynamics of the 



15 



equations when (3 = is reflected by the fact that the splay state has N — 2 unit multipliers 
when (3 = 0. Conversely, if the dynamics is nondegenerate, then there must be exactly one 
unit Floquet multiplier. 

We calculate the multipliers of the splay state as a function of (3 for two types of the 
shunt loads, an LC load and a purely capacitive load. We find that there is a single unit 
Floquet multipler when (3 > 0, indicating that the system is nondegenerate. However, when 
(3 is larger than about 1, there are N — 2 multipliers very close to 1. This is why we observed 
an apparently neutrally stable incoherent state in Fig. 0(d): the Floquet multiplier which 
measures the speed of the drift in the incoherent manifold is extremely close to 1. 

Furthermore, a complex conjugate pair of multipliers approaches 1 as (3 — > oo. Thus, in 
the large f3 limit there are N unit Floquet multiplers, corresponding to the oscillators being 
effectively uncoupled, which we observed in the numerical integrations. 

In order to calculate the Floquet multipliers, we first find a stable splay state for one set 
of parameters, then use a continuation method to compute the branch of periodic solutions 
as (3 is varied, holding the other parameters fixed. Solutions can be obtained regardless of 
stability, which clarifies the dependence of the multipliers on (3. 

The starting solution is found by choosing the same parameters as two of the cases studied 
by Nichols and Wiesenfeld [NicfcWie92|| , specifically, their cases II and III: 

[NWII] N = A L = 0.75 R = C = 20 I b = 2.5 (3 = 1.1 
[NWIII] N = 4 L = R = C = 1 I b = 1.5 (3 = 0.1 

(We have scaled their actual parameter values to suit our definition.) 

In both cases, we have chosen an initial condition close to the splay state. Since the 
waveform is nonuniform, <f)j(0) = 2nj/N is not a good starting point. Instead, we have 
computed the periodic solution <f>*(t) of the uncoupled equation 

f3<j)j + 4>j + sin <fij = I b (5) 

for the given (1^,(3), then assigned the initial condition as 0^(0) = <p*(Tj/N) where T is 
the period. The choice is especially useful for the LC-load (case II), and it enables us to 
obtain the splay state (or, at worst, an incoherent state extremely close to the splay state) 
sufficiently fast. We do not need an extrapolation technique as used in ||NicfcWic92|| . 

Next, we sweep (3 up to 10 using the continuation and bifurcation package AUT086 
PoefcKer86|| . Then, (3 is decreased from 10 down to almost 0. AUT086 computes the 



Floquet multipliers, but not accurately enough for our purposes. To improve the accuracy, 
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Figure 6: The magnitude of Floquet multipliers as a function of (3, for (a) an LC load (case 
[NWII]) and for (b) a C load (case [NWIII]). All multipliers are shown in the left figures while 
the shaded regions are enlarged in the right figures. In order to resolve the real multiplier 
near unity, another set of enlargements (of the shaded region in the right figures) are shown 
as insets. There are 10 multipliers for (a) and 8 for (b). The letters denote "cc" (complex 
conjugate pairs), "r" (real), and "u" (the unit multiplier). The arrows labeled "NW" point 
to the (3 values used in |[NicfeWie92l . 



we output a point on the periodic solution at each value of (3, and integrate the Jacobian 
matrix for one full period. The eigenvalues of the resulting matrix are the Floquet multipliers, 
which we compute using Mathematica. There must always be a unit multiplier corresponding 
to perturbation in the direction of the periodic orbit. This multiplier is computed to 8-12 
significant figures. (AUT086, which needs Floquet multipliers only for detecting bifurcation 
points, has calculated the trivial multiplier within 3-7 digits. An improved version AUT094 
was not available at the time of our computation.) Also, the multipliers for the starting 
value of j3 agree with those obtained in [ [NicfeWic92|| . 

Figure |6] shows the magnitude of the computed multipliers. The figures (a) and (b) 
correspond to [NWII] and [NWIII], respectively. The left figures show all the multipliers 
whereas the right ones are enlarged views of the shaded box in the left figure. To resolve 
the real multiplier near unity, the shaded box in the right figure is further enlarged as the 
inset. The number of multipliers is the same as the dimension of the system, thus there are 
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ten multipliers in (a) and eight in (b). They are labeled as "cc" (complex conjugate pairs), 
"r" (real), and "u" (the unit multiplier). One sees that the unit multiplier is computed 
accurately, and shows no visible deviation from unity even in the inset. 

Although these multipliers provide only local information around the splay state, we can 
roughly associate them with the global dynamics described in Section 2. For example, five 
multipliers in (a) and four in (b) have magnitudes well below unity, especially when (3 is 
small. They correspond to the rapid initial transient when the junctions shuffle and change 
firing order. After these directions are contracted, all the junctions have similar periods. The 
convergence speed of this initial stage is instantaneous when (5 — and becomes slower as f3 
becomes large. One real multiplier in (a) which always stays near 0.9 seems to correspond 
to a perturbation in the charge Q on the capacitor. This multiplier behaves differently from 
the others. 

After the initial transients are over, the time-delayed coordinates described in Section 2 
become meaningful, and the dynamics involves the phase differences. The three nontrivial 
multipliers of near- unit magnitude, shown in the right panels in Figs. [](a,b), represent this 
dynamics. From now on, our interest is focused on these three multipliers. The incoherent 
manifold is the key to understanding the dynamics. Two multipliers (the cc. pair) corre- 
spond to the dynamics transversal to this manifold. If they have magnitude greater than 
unity, then trajectories depart from the manifold and the in-phase state is usually stable. 
When their magnitude is below unity, the incoherent manifold is attracting, at least near 
the splay solution. 

The real multiplier near unity, which we call the "critical multiplier," represents the 
dynamics along the manifold. (We describe 4 oscillators. In general there are N — 3 critical 
multipliers.) When the critical multiplier is less than 1, the splay state is attracting within 
the incoherent manifold. When it is larger than 1, the drift on the incoherent manifold is 
away from the splay state. In the (3 = equations the critical multiplier is exactly unity, 
and there is no drift on the incoherent manifold. 

Figure |B| shows that the magnitude of the "transverse multipliers" deviates substantially 
from 1 when f3 is small. The magnitude changes in a nontrivial fashion, switching stability 
at finite f3 in either case. Using an LC load, the incoherent manifold is repelling for small 
j3, and attracting for larger values. For a C load, it is reversed. 

The critical multiplier is extremely close to (but not equal to) 1 at moderate j3. Thus it 
appears that degenerate phase space structure persists for j3 > ||NicfeWie92|| , since the drift 
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on the incoherent manifold is too slow to observe. However, we believe that the dynamics 
is nondegenerate for (3 > 0. At those (3 values where the critical multiplier crosses unity, 
we expect a pitchfork bifurcation that creates two generic incoherent solutions. Analysis 
WatfeStr94ll predicts that the real multiplier must be exactly 1 at (3 = 0. Fig. | is in 



agreement. (Yet another enlargement of the inset of Fig. 6(a) shows that the critical 
multiplier dips below 1 at very small f3 and then approaches 1 as (3 — > 0.) 

4 Analysis 

In this section we show analytically the weak destruction of the foliated phase space described 
in Sec. 2, and estimate the Floquet multipliers around the splay state computed in Sec. 3. 
The analysis uses the multiple time-scale method, and is applicable for a general N, but our 
primary focus is on the case N = 4 in order to compare with the previous sections. 

We introduce a small parameter e and different time scales in Sec. 4.1, then study the 0(1) 
equations in Sec. 4.2. This lowest order system describes the dynamics when the junction 
phases are being shuffled and the time-delayed coordinate is not yet meaningful. In Sec. 4.3 
we go up to the 0(e 2 ) expansion, and derive an averaged equation, which describes the drift 
of the phases. The order parameter Ri is shown to converge toward either or 1. The 
averaged equation can predict the magnitude of the transverse multipliers, but the critical 
multiplier is not yet resolved. Thus, in Sec. 4.4 we go up to 0(e 4 ), and derive a higher order 
averaged equation, which is capable of describing the drift along the incoherent manifold. 
Now, i?2 converges toward or 1. This predicts the last remaining multiplier when N = 4. 
In Sec. 4.5 we summarize the results. Since the calculations are standard, details are shown 
in the Appendices. 

4.1 Perturbation expansion 

Chernikov and Schmidt | ChcfcSch95 l recently used the first-order averaging method to study 



the governing equations Eq.(|J) with f3 > 0, and obtained the condition for the junction phases 
to synchronize to the in-phase state. It is not clear how to extend their method to higher 
order. Here, we employ a multiple time-scale analysis that allows us to go up to higher order, 
and extracts more information on the asymptotic dynamics and the phase space structure. 
We first rescale the time from t to r as 

t = r/I b , (6) 



19 



and rewrite the governing equations Eq.(|I]) as 



rri(j)j" + (f)j + e sin <j>j + Q' = 1 , j = 1, . . . , N, 

1 N (7) 
LI b Q" + RQ' + Q/{CI b ) = - Y, <t>k { } 

iV k=l 

where m = 1^(3, e = 1/h, and primes denote derivatives with respect to r. Following 
|ChefcSch95|1 , we fix m, and study the limit e — > 0+. In the pendulum analog of the 



Josephson junctions, m is a scaled mass of the pendulum and e measures the ratio of the 
gravity force to the applied torque. 

As usual, the perturbation expansion will generate secular terms as we proceed to higher 
order. In order to eliminate them, it turns out that we need to introduce following three 
different time scales: 

T = r, T 2 = e 2 r, T A = e\ (8) 
which will be formally treated as independent variables. Then, 

' = d + e 2 d 2 + e 4 d 4 

" = d 2 + 2e 2 d d 2 + e 4 {d 2 2 +2d d i ) + O{e 6 ) (9) 
where d = d/dT etc. We expand the solution in powers of e: 

oo oo 

^(r;e) = £ £-^(T , T 2 , T 4 ) , Q(r;e) = £ e p Q p (T , T 2 , T 4 ), (10) 

p=0 p=0 

as well as the sinusoidal nonlinearity: 



sin < 



3 



Ys p S pj (11) 

p=0 

where S 0j = sin0 O j, Sij = y cos0 O j, S 2j = 2 jCos0 O j - (<t>lj/2) sin <f) 0j , S 3j = 4> 3j cos (fi j - 
(0ij/6) cos^oj — 4>ij<f>2j sin0 O j, and so on. All (p P j and Q p are assumed to be bounded except 
4> j, which increases without bound in time. Without loss of generality, we impose that the 
averages of <p P j, p > 0, over the fast scale T vanish. 

4.2 Fast dynamics 

Substitute Eqs.(p|-pJ]) ©> and equate like powers of e. At the lowest order 0(1) we 
obtain 

mdl(j) Qj + d 4> 0j + d Q = 1 , j = 1, . . . , N, 

1 N (1?) 
LI b d 2 Qo + Rd Q + Q /(CI b ) = -Y, d ° ( l ) °k. 1 ; 

k=l 
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(a) LC-load 



(b) C-load 




Figure 7: For (a) LC-load and (b) C-load, the data points in Fig. 6 whose magnitudes 
are well below unity are estimated by the 0(1) expansion. Solid curves are magnitudes 



corresponding to Eqs. ([l4| , |T5D . 



This system is linear, inhomogeneous. Its periodic solutions are easy to find. 



T + 8j , Qq — Ch 



(13) 



where 9j are arbitrary "constant" phases which may, in fact, depend on T 2 and T 4 . When all 
Oj are identical (or equally spread apart), the solution is the in-phase (or splay) state. The 
(residual) phases 9j are same as the phases obtained by the time-delayed method of Sec. 2. 

As shown in Appendix A, it is easy to see that any initial condition converges to a solution 
of the form ([13]), according to the flow flT2|). The convergence is exponentially fast in the 
time scale of To, and the rate predicts some Floquet exponents. In particular, the periodic 
solutions ([□]) must have (regardless of the choice of 9j) real negative Floquet exponents 



A 



-27c/m (with multiplicity N — 1). 



Also, another set of exponents are obtained by solving 



mLI b X 3 + (mR + LI b )\ 2 + 



f m 



+ R+1 A 



0. 



(14) 



(15) 



\CI b CI b 
This characteristic polynomial produces 1-3 roots depending on the type of the load, but 
their real parts are always negative. 

Thus, all the exponents except N are estimated at this order. In Fig. 7 they are converted 
to multipliers \i = exp(27rA) (since the period is 2%), and shown as solid curves. (The Floquet 
multipliers are invariant under the scaling of time in (0).) The agreement with the data points 
located well below the unit magnitude is excellent even though the formal perturbation 
parameter e = 2/3 is not so small. These rapidly converging multipliers correspond to the 
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initial shuffling of oscillators, and after these directions are contracted, the solution is of the 
form ( P~5| ) so that all the junctions have the same period. 

The remaining N multipliers correspond to perturbation in the phases 8j (j = 1, . . . , N). 
At this order of expansion the junctions are frequency locked, but there is no preferred 
phase, so that any perturbations on 9j remain as applied. Therefore, the N multipliers 
are all estimated to be unity, and we need to go up to higher order to resolve the issue of 
neutrality. 

4.3 Transverse multipliers 

The arbitrary phases dj are actually dependent on T2 and T4. The drift of the phases can 
be studied by assuming that the fast dynamics has already contracted, and <f) j and Qq are 
given by (|IBD. Geometrically, we are only interested in the dynamics on the center manifold 
(the invariant iV-torus) of the system (|12|). 

Since we find no secular terms at the 0(e) expansion, we need to consider up to 0(e 2 ). 
(See Appendix A.) There, we must impose a solvability condition to prevent 4>2j from growing 
unbounded. The condition is found to be 

1 N 

d 2 9j = ^ + - £ K cos(0 fc - 9 3 ) + h sm(6 k - 6,)} (16) 
fc=i 



where 



and 



Ql = o7TT~ 2^ > fl i = 2Re N > 6 i = -2Im[ Cl ] , Cl = 1 (17) 

2(1 + m z ) 4(1 + im)Xi 

Xi=(l + im)Z 1 + 1 , Z± = R + i ( LI b ' 



CI b 

Note that Zi is the impedance of the load for the frequency lb- We have not found a physical 
interpretation for X\ or c\. 

Equation (|16D describes the slow drift of the phases. It was derived by [|ChefcSch95 



using the averaging method. It gives the 0(e 2 ) contribution to the normal form for the 
dynamics on the center manifold, which is the invariant iV-torus. We call this normal 
form the averaged equation, since at each order we remove the next harmonic of the phase 
modulation superimposed on the fundamental oscillation. See |[AshfcSwi92| | . Thus, the 
second order averaged equation is 



63 = 1 + e 2 (fti + ^ E N cos (^ " + b i sin (^ - 0j)>) + 



(19) 
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where the dot represents d/dr. This has the same form as the averaged equation obtained 
in the overdamped limit (3 — > of the governing equation ([I]) | |Jaifeal84| , |WiefeSwi95| . 



Strictly, the second order equation for 9j would look like Eq. (|T^), except without the 
"1", which represents a constant frequency. The equation for O j = T Q + 8j has exactly the 
form of Eq. ([19]). For simplicity, and agreement with the notation in ||AshfeSwi92|| , we have 
replaced (fioj with 9j in Eq. ([OJ). 

Equation (|lj]) is fully nonlinear, but is solvable by use of a coordinate transformation 
WatfeStr94|| . It is proven that the sign of the parameter b\ determines the attracting set; 
when bi > the phases 8j eventually converges to the in-phase state while for b\ < 
the phases tend to spread so that they converge to one of the incoherent states. The first 
moment R\ defined in Eq.(|3D, monotonically converges either to R\ — 1 (in-phase) or R\ = 
(incoherent), depending on b\. The critical case b\ = has been termed by ||ChefeSch95|| the 
"synchronization condition" . In our notation it reduces to 

(1 - m 2 ) (LI b - — ) + m{l + 2R) = 0. (20) 

Floquet multipliers of the system ( |T9"D can be evaluated easily around both the in-phase 
and the incoherent states. (See ||AshfcSwi92|| and Appendix B here.) The in-phase state has 
the unit multiplier, along with a real Floquet multiplier with multiplicity N — 1: 

^- phMe = eXP ll + ^(fi 1 + a 1 )J (21) 

Assuming iV > 3, the splay state (6j evenly distributed) in equation fll9|) has a unit multipier 
with multiplicity N — 2, along with a complex conjugate pair of Floquet multipliers with 
magnitude 

M = ap {TT^o:)- (22) 

These approximate the transverse multipliers we computed numerically in Sec. 3. Fig. 8 
shows a comparison of Eq.(|22"D (solid curves labeled as "T") with the numerical data from 
Fig. 6. There is a systematic deviation in the magnitude, probably due to a large value 
of e. Nonetheless, we see that the formula (|2"2] ) provides a good estimate of the transverse 
multipliers. In particular, the transition value of (3 when |/i| crosses unity is estimated 
accurately. 

When bi > 0, the analysis concludes that the in-phase state is the attractor, and higher 
order corrections will provide little information. However, when b\ < 0, the truncation in 



23 



1.1 



1.0 



(a) LC-load 




1 003 












| 1.000 

\ / 


/ \ 0.8 


1.2 


0.4 


(3 

P 




4 




Figure 8: Comparison of the data points in Fig. 6 with the perturbation analysis. The 
curves labeled as "T" are the estimate (B^) of the magnitude of the transverse multipliers. 
The curves labeled as "C" in the insets are the formula ( p8|) for the critical multiplier. 

equation (plj) is too severe. The splay solution is stable in the transverse direction, but there 
are N — 2 unit multipliers whereas there is generically a single unit multiplier. The truncated 
system has an attracting (N — 2) -dimensional torus, the incoherent manifold, foliated with 
periodic orbits. Higher order terms break this foliation and lead to generic dynamics. 



4.4 The critical multiplier 

To eliminate the unit multiplers caused by the truncation of equation (JIP|), we must go to 
the next order. By imposing a new solvability condition, we eliminate the secular terms 
which arise at the 0(e 4 ) expansion (since there are no such terms at 0(e 3 )). The calculation 
is shown in Appendix A. We find the solvability condition to be: 

1 N 



d^j =n 2 + -Y,{a 2 cos 2{6 k - 9j) + b 2 sin 2(9 k - 6 3 )} . 
^ k=i 



Here, 



15m 4 + 10m 2 + 1 
'8(1 + m 2 ) 3 (l + 4m 2 ) 



and a>2 = 2Re[c2], 62 = — 2Im[c2] where 

c 2 

and 



vm 



16X2(1 - im)(l + im) 2 (l + 2im) 



X 2 = (l + 2im)Z 2 + 1 , Z 2 = R + i[2LI b - 



2CI h 



(23) 



(24) 



(25) 



(26) 
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Note that Z 2 is the impedance of the load at frequency 2I b . 

The drift occurs at the time scale of 0(e 4 ), and the fourth order averaged equation is 



As described in Appendix B, the order e 4 terms we calculated contribute to the stability of 
the splay state. If N > 4 there is a complex conjugate pair of Floquet exponents proportional 
to c 2 and ci. For iV = 4, the single critical Floquet multiplier is 



In the insets of Fig. 8 the predicted critical multiplier is plotted for both types of the load 
(curves labeled as "C"), and comparison is made with the numerical data from the insets 
of Fig. 6. Again, we observe a systematic deviation, but dependence of the multiplier on f3 
is qualitatively well predicted. The critical (3 where the splay state changes its stability is 
estimated accurately. 

4.5 Global dynamics and attractors 

The fourth order perturbation expansion allows us to predict the attractor for the Josephson 
junction system with N = 4 oscillators. The incoherent manifold (or incoherent bar) defined 
by Ri = is dynamically invariant when N = 4 (see Appendix B), and we can discuss the 
drift on the incoherent manifold. For N > 5, the set R± = is not dynamically invariant. 

We have concentrated on the Floquet multipliers of the splay state, but the stability of 
the in-phase and the 2+2 incoherent states can also be calculated. We find that, for N = 4 
and sufficiently small e: 



This table does not show the bistability or bi-instability of states, which depends on 
corrections to the Floquet exponents at higher order in e. For example, when the order e A 
terms in the averaged equations are truncated, the Floquet exponents of the splay and 2+2 
incoherent states are the same, and these have opposite sign of the Floquet exponents of the 




(27) 



+ O^e 4 ) + 0(e 6 ) 




(28) 



If bi > then the in-phase solution is asymptotically stable 

If b\ < and b 2 < then the splay solution is asymptotically stable 

If b\ < and 6 2 > then the 2+2 i solution is asymptotically stable 



(29) 



25 




1 2 3 I b 2 4 I b 



Figure 9: The asymptotically stable state predicted for N — 4. The perturbation theory 
assumes e = < 1, In Fig. (a), the same load as in Fig. 2 (L = 1, R — 0, C — 1/4) is 
used. Fig. (b) has L = 100/9, i? = 0, and C = 1/100. The dots in Fig. (a) marked as "a" 
through "d" correspond to the parameter values used in Fig. 2, respectively. The solid curve 
is where bi = 0, and the in-phase state is stable in the shaded region. Outside the shaded 
region, the incoherent manifold is attracting. There is a slow drift on the manifold whose 
direction switches on the dashed curve 62 = 0. The splay solution is asymptotically stable 
in the upper-right region, but the drift is extremely slow. The point "b" where we observed 
a stable splay state in Fig. 2 is located in the region for the 2+2 incoherent ("2+2 i") state; 
we attribute this to the rather large value of e. 

in-phase solution. All of these exponents are proportional to 61. However, when order e A 
terms are included in the averaged equations, the exponents for each of the three solutions 
gets a different correction. 

Similarly, when the order e 6 terms in the averaged equations are truncated, the drift on 
the incoherent manifold is toward the splay solution (R 2 — > 0) if b 2 < and the drift is 
toward the 2+2 incoherent solution (R 2 — > 1) if b 2 > 0. However, when order e 6 terms are 
included, the exchange of stability is through a generic incoherent state that is born and dies 
at separate pitchfork bifurcations of the splay and 2+2 incoherent states. 

To summarize, the asymptotic analysis shows that R\ and R 2 monotonically changes in 
separate time scales in the limit e — > 0. For N = A this determines the attracting periodic 
solution uniquely for a given set of parameters; the possible ones are the in-phase, the splay, 
or the 2+2 incoherent state. Fig. 9 shows the attractor we predict in the (h, P) plane, using 
two different LC loads. 

The synchronization condition b± — determines when the in-phase state is stable. By 
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substituting m = hfi into (p0[) , we obtain a quadratic equation for I 2 : 

(LCI 2 - 1)((3 2 I 2 - 1) = /3C(1 + 2/2) (30) 

For L > 0, there are (at most) two positive critical values for each f3. The boundary is 
shown as solid curves, and the region where the in-phase solution is stable is shaded. When 
P = 0, b x < for I b < (LC)- 1/2 and bx > for I b > (LC)- 1/2 . This agrees with ||WiefeSwi95 



and allows us to determine the sign of b\ in each region. When C is small, the two factors 
on the left-hand side of Eq. ( p0|) determine the asymptotic boundaries seen in Fig. 9(b). 

The dashed curve, where 62 = 0, concerns the direction of flow in the incoherent bar. 
The condition also becomes a quadratic equation for I 2 : 

(ALCI 2 - l)(4/? 2 / 6 2 - 5) = C[4/3(l + AR) - 2(1 + R)/{3] (31) 

Furthermore, 62 = when /3 — 0. Sign of 62 can be determined from the fact that 62 < for 
P sufficiently small and positive. When C is small, the two factors on the left-hand side of 
Eq. ( |31"D determine the asymptotes to b 2 = seen in Fig. 9(b). 

A plot similar to Fig. 9 can be made for any load (not just an LRC load) by using the 
impedance of the load at frequency and 2I b for Z x and Z 2 in equations ( |18D and (|26|). 

The bistability we observed numerically in Sec. 2 is not computed by our analysis. As 
described after Eq. (p9[) , higher order corrections split the curve solid curve ib\ = 0) into 
three curves, and the dashed curve (62 = 0) splits in two. 

In Fig. 2(d) we used = 2.5 and /3 = 1 and observed that the system apparently con- 
verged to a generic incoherent state. The splay is asymptotically stable here, but the drift 
toward it is too slow to observe. Farther into the upper-right "splay" region of Fig. 9(a), the 
coupling between the oscillators is so weak that even the dynamics transverse to the inco- 
herent manifold becomes too slow to observe. The oscillators "fall" into a phase determined 
by the initial conditions, and then Ri and R2 appear to remain constant. 

The lower-left region where the in-phase is unstable is the most dynamically interesting, 
since this is where the moments R\ and R2 evolve quickly. In the region near I b = 1 there 
are period doubling bifurcations and bifurcations to fixed-points and semirotors |[Arofcal95 



which are not shown in Fig. 9. The asymptotically stable splay solutions and 2+2 incoherent 
solutions exist in the regions with I b > 1. 

In Sec. 2.2 we describe the numerical observations along the line l b = 1.5. There is 
qualitative agreement, but the transition between splay and 2+2 incoherent occurs at /3 ~ 
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0.2, higher than in Fig. 9(a). This discrepancy is not surprising, considering that h = 1.5 
is not large. (On the other hand, recall that the (3 value where h 2 = is predicted quite 
accurately in Fig. 8, where h = 2.5 and = 1.5 with different loads.) The solid curve in 
Fig. 9(a) crosses lb = 1.5 within the observed range of bistability between 2+2 incoherent 
and in-phase (0.40 < (3 < 0.61). 

We have compared our prediction with the numerical computations of Hadley et al. 
[ Hadfcal88|l for the stability of the in phase state. The agreement is pretty good, especially 



for their Fig. 2(e), where L = C = 1/2, and R = 0. Our figure for this load looks similar 
to our Fig. 9(a), since the resonant frequency {LC)~ 1 ^ 2 = 2 is the same for both loads. The 
stability boundary of the in-phase that starts at lb = 2 and arches to the left of lb = 1 agrees 
surprisingly well with [ Hadfcal88| . (This is why we show the region with lb < 1 in Fig. 9.) 



For the other LRC loads considered in ||Hadfcal88| , our agreement was good for lb larger 
than 1 or 2, but not good for smaller I b . 

5 Summary and discussion 

Globally coupled oscillator systems arise not only in the Josephson junction context but also 
in lasers |0ts91| , [5il&al92 \ Rap_93], the complex Ginzburg-Landau equation | HakfeRap92 



[NakfcKur94|| , and coupled maps |[Kan91|| . A ubiquitous feature among these systems are the 
apparent neutral stability of the splay state over a wide range of parameters. In a model 
of a laser array [ Rap93 ] it is trivial to see that a foliated incoherent manifold exists, and 



hence that the splay solution is neutrally stable. Series Josephson arrays with f3 = have a 
similar structure, but it is more difficult to prove ||WatfeStr94]| . On the other hand, for series 



Josephson arrays with (3 > we have shown that the incoherent manifold is not foliated by 
incoherent periodic solutions, and that splay solution is not neutrally stable. 

We have also explained why previous numerical studies appear to indicate that the splay 
solution is neutrally stable even when (3 becomes finite: The approximate neutral stability 
is seen when (3 and/or lb are large. When lb and (3 are moderate, we can easily find asymp- 
totically stable splay solutions for an LC or C-load. We have visualized the global dynamics 
and the drift on the incoherent manifold (bar) when N = 4. Our multiple-scale analysis 
described this process quantitatively. 

Recall that e — 1/ 'lb and m = Ib(3, and our perturbation expansion assumed m was fixed 
while e — *• 0. This is the natural expansion given the form of the equations, but a different 
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scaling of the equations shows that another weak coupling limit is (3 — > oo with fixed. 
On the other hand, (3 — > with If, fixed is not a weak coupling limit, because £1% and c\ in 
Eq.(17) have a nonzero limit as m — > 0. 

The previous discussion of weak coupling limits makes no assumptions about the load. 
Another way to get weak coupling is to make the load impedance large. For example, 
R — > oo is a weak coupling limit, as exploited in ||WiefcSwi95|| using a different scaling of the 
Josephson junction equations. 

The approximate degeneracy which occurs in series Josephson junction array is due to the 
simplicity of the averaged equations when truncated to order e 2 . For N oscillators, we must 
go to order N (if N is even) or iV — 1 (if iV is odd) to remove the degeneracy of the truncated 
averaged equations (see Appendix B). For N = 4 junctions, the 0(e 2 ) equation (|T6|) was not 
sufficient, but adding the 0{e A ) corrections in equation (|23|) remove the degeneracy. 

Apart from the theoretical interest on whether the neutrality is exact or not, our anal- 
ysis re-confirms the difficulty of using the splay state for practical purposes. Firstly, the 
parameters must be chosen so that the state is attracting. For N = 4 the two conditions 
b\ < and 62 < are required; for a larger N, there will be more conditions to satisfy, 
which would make it necessary to tune the parameters precisely. Secondly, even when this 
is possible, the convergence would be very weak. For N = 4, the drift of R2 is slow unless 
the capacitive load is used and the parameter values are chosen carefully. For a larger N, 
the drift of the higher phase moments R n occurs at a rate proportional to e 2n , and the splay 
state is practically neutral in many directions when I^ = 1/e is large. 

It is unknown what happens for large N when 1^ is near 1. For example, if the parameters 
are chosen as in Fig. 2(b) but N is large, is the splay still stable? This depends on the signs 
of b n , for n up to floor [iV/2]. Based on how b\ and 62 depend on the impedance of the load 
at frequency I b and 2I b , respectively, we expect that b n is positive but near for large n, 
since the impedance of the load at frequency nib is large and inductive. This will cause the 
splay to be unstable, with approximate neutral stability. On the other hand, it is possible 
that a capacitive load has a robustly stable splay solution for any number N of Josephson 
junctions. 

In order to get a robustly stable splay solution for practical applications, (3 should be 
small but nonzero, and h slightly larger than one. Furthermore, a C-load should be used 
(in any case L should be a few times smaller than 1/C, although the coupling will be weak 
if 1/C is too large). In a different approach to make the splay state more stable, recent work 
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Dar&al96 



|Har&Gar96| 



describes modifications of the series array circuit. 
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A Perturbation expansion 

In this Appendix we show the detail of the calculation in Sec. 4. 
A.l The fast dynamics 

We start from the setup described in Sec. 4.1. The 0(1) equations are Eq. (|T2|) . The stability 
of the periodic solutions (|T3"D can be studied by plugging 



? 0j = <P Qj 



+ Uj , Q = Ql + q, (32) 



the choice of 9j, we obtain 



into (|12|) and writing an equivalent linear homogeneous system for Uj and q. Regardless of 

mdluj + d Uj + d q = , j = 1, . . . , N, 

LI b dlq + Rd q + q/(CI b ) = ^ £ d Q u k ' (33) 

k=l 

There are as many independent solutions as the system dimension (2iV + 2 for Li?C-load 
and 2N for C-load) that grow by a scalar factor \x = exp(27rA) after one period 2ir in T . 
The growth factor /j, (or A) is a Floquet multiplier (or exponent). Since the system is time 
independent, A is simply an eigenvalue of ([!3|) and is easy to find. 

First, there are iV solutions corresponding to perturbations in phases only: 

Uj = 5 jn , d Uj = , q = d q = (34) 

for n — 1, . . . , Af (where 5j n is the Kronecker delta function). Each n corresponds to one 
unit multiplier \i = 1, A = 0. 
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The second group corresponds to the solutions of the form: 

Uj = Aj exp(ATo) , q = , A = — 1/m (35) 
which satisfies the first equation in Eq.(|33|). The second equation can be satisfied by imposing 

N 

£ A k = 0. (36) 
fe=i 

Since there are N parameters Aj with one constraint fl36|), there are N — 1 Floquet multipliers 
/x = exp(— 2ix/m) from this group. 

The last group of multipliers come from the solutions of the form: 

Uj = Aexp(XTo) , q = —A(mX + 1) exp(AT ) (37) 

where A is non-zero. Again, the first equation in Eq . (|33"D is satisfied, and the second equation 
can be also satisfied by imposing Eq.fll5|). The characteristic polynomial for A has 3, 2, and 
1 root(s) for an LRC-, RC-, and C-loads, respectively. Using Hurwitz's criterion, it is easy 
to show that the real parts of the roots are always negative. 

A. 2 The first solvability condition 

At 0(e), we obtain 

mdlfaj + dofaj + d Qi = - sin Oj , j = 1, . . . , N, 



LI h dlQ x + RdoQi + Qi/(CI b ) = i- d <p lk . ^ (38) 

iV k=i 

The system is linear, with a forcing from 0oj- The homogeneous part for <f>ij and Q\ is 
identical to that of the 0(1) expansion for O j an d Qo- Thus, the homogeneous solutions 
consist of exponentially decaying terms with possible constant phases to <p\j. We set the 
constants to zero so that the average of <p\j over T becomes zero. 
At this order no secularity arises. The solution is 



> Xj = y A X j (T 2 , T 4 )e iTo + c.c.) /2 , j = 1, . . . , N, 
Q 1 = (B 1 (T 2 ,T 4 )e iTo +c.c.)/2 



(39) 



where "c.c." stands for the complex conjugate, and 



1 / 1 v 

1 I AO. 1 



i N 

1 V AT ^ 



(40) 



X * N k=i 
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with X\ given in (|I8"D. 
At 0(e 2 ) we obtain 

mdlfaj + do<p 2 j + d Q 2 = -d 2 (p j ~ <Pij cos</> 0i , j — 1, 



iV, 



LI b dlQ 2 + i2a Q 2 + Q 2 /(C/ 6 ) = jj E(fib^ + 9 2 ofc ) 

fe=i 



(41) 



The forcing term on the right hand side of the first equation involves "dc" components 
(constant of T ). Unless their sum vanishes, 2j - becomes unbounded in T . Thus, we impose 
the solvability condition flTBf). 

The solution is then (neglecting the homogeneous part as in the previous order) 

2j = (A 2j (T 2 , T 4 )e 2tTo + c.c.) /2 , j = l,...,N, 
Q 2 = Q 2 o(T 2 ,T 4 ) + (B 2 (T 2 ,T 4 )e 2lTo + c.c.) /2 



(42) 



where 



.4 



2.7 



z/4 



A' 



Q 



1 + 2im 



A Xj e 1 ^ 



CL N 



J = l 



,...,iV, 



20 



X (A lfc e-** + c.c. 



Bo 



(43) 



with X 2 given in (^g 



A. 3 The second solvability condition 

We need the higher order expansions only to resolve the dynamics along the incoherent 
manifold. On this manifold we have 



N 



k=l 

identically. Thus, we simplify the terms accordingly. 



(44) 



A Xj = e^/(l + im) , B x = , Q x = 
= -1/2(1 + m 2 ) = , Q 20 = ChQ, 



A 2i 



i/A 



f l+im)(l + 2im) 
i/A 



e 2idj _ 



X 2 N 



E 



Bo 



-Ye 2 
1 + im)X 2 N ^ 



(45) 

(46) 
(47) 

(48) 
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Now, we write down the 0(e 3 ) expansion. 

' mdftfoj + dofoj + d Q 3 

= -2md d 2 (t)ij - d 2 (f>ij - <p 2j cos (f> 0j + {<j>\J2) sin (p 0j , j = 1, . . . , N, 

i n (49) 

LI b d 2 Q 3 + Rd Q 3 + Q 3 /(CI b ) = -]T 9 o 3fc 

fc=i 

The forcing terms in the right hand side only involve exp(±zTo) and exp(±3zTo) so that 
there is no secularity. The solution is 



(A 31J (T 2 , T A )e iTo + A 33J (T 2 , T 4 )e 3iTo + c.c.) /2 , j — 1, . . . , N, 
Q 3 = (B 31 (T 2 , T 4 )e lTo + B 33 (T 2 , T A )e 3tT ° + c.c.) /2 



(50) 



with 



■n 2 + llim — 11m 2 

A 31j e~^ = 



[1 + im) 3 (l — im)(l + 2im) 



' j^e 2 ^-^ (51) 



5(1 + im) 2 (l + 2im)X 2 N ^ 



and 



B 31 = 0. (52) 

(We do not need to write down A 3S j and -B33.) 

Finally, we consider the 0(e 4 ) expansion. Here, terms independent of To appear in the 
equations for faj, and we will set these terms to vanish. This condition becomes 

d A 9 j = - (A 31j e~ i(> i /4 - .^ ~\ : ,< '" /32 - /7T| r l i; < '"• /8) + c.c. (53) 

which reduces to Eq. (|23|) . 



B Averaged equations and Floquet exponents 

This Appendix concerns the structure of the averaged equations, and shows how to pick out 
the terms that determine the stability of the splay solution. 

We also describe three special types of averaged equations: pairwise coupled systems, 
centroid coupled systems, and "the" averaged equation which is both pairwise coupled and 
centroid coupled. Contrary to popular belief, none of these special systems describe the 
weak coupling limit of the series Josephson junction arrays with (3 > 0. We conjecture that 
a centroid coupled system is obtained in such a limit for the array with (3 = 0. 
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We have observed that the perturbation expansion described in Section 4 and Appendix 
A leads to an averaged equation of the form 

9j = l + e 2 [a + (ae~ i ^z 1 + c.c.) 

b + c\z x \ 2 + (e- ie i((3zi + 7^iN 2 ) + e" 2 ^(/iZi 2 + uz 2 ) + c.c.)] (54) 



+ e 4 



+ 0{e b ) 

where a, b, and c are real constants and a, /3, etc. are complex constants. The permutation 
invariant z m are the moments of the N oscillators placed on the unit circle, defined in 
equation (|3|). 

We conjecture that if the the calculation were taken to arbitrary order, the resulting 
averaged equation would have the form: 

9j = 1 + F(ee ie i,ez 1 ,e 2 z 2 ,e 3 z 3 ,...), (55) 

where the function F in equation ( [55] ) is real- valued and invariant under the phase shift: 

e .^ ej + t, z m ^z m e imt (56) 

See section 2 of | AshfcSwi92 | for a discussion of this phase shift symmetry, which is analogous 
to the phase shift symmetry of a Hopf bifurcation normal form. Properties of the dynamics 
that depend on the phase shift symmetry do not carry over exactly to the unaveraged 
equations. Furthermore, F has the property that each monomial term has at least one 
factor of ee ±j£lj , possibly in the form (ee j )(ee~ ie: >) = e 2 . For example, the term £ 2 |zi| 2 is 
found on the right-hand side of Eq. fl5"4|). The lowest order term of this form is e 4 |2;i| 2 , with 
coefficient c. 

We cannot prove that the perturbation expansion to any order leads to a system of the 
form (|55|) . The reason for our conjecture is this: For the Josephson junction system, the 
coupling of the oscillators is through the current Q, which is permutation invariant and can 
be written in terms of the z m . There must be at least one exponential factor (e ±l ^) for the 
coupling back to oscillator j. The ordering of the terms as powers of e follows because each 
exponential factor e l9k is tied to a power of e. For example, z 2 is quadratic in e* £lfc , and it 
appears with a factor of e 2 . We expect this to continue to higher orders due to the sin(</>j) 
type nonlinearity of equation ([!]). 

The truncation of terms of order e 4 and higher in ([54]) leads to a system so common that 



we have nicknamed it "the" averaged equation. See ([19|) for "the" averaged equation of the 
Josephson junction system. 
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The dynamics of "the" averaged equation are well understood, and highly degenerate. 
First of all, there is an invariant incoherent manifold foliated with periodic orbits. Thus, the 
splay solution is neutrally stable to all perturbations within the invariant manifold. Secondly, 
the stability of the in-phase solution and the normal stability of the incoherent manifold are 
both determined by the sign of b\. If one solutions is stable, the other is unstable: there can 
be no bistability. As a result, "the" averaged equation does not have a generic bifurcation 
as bi passes through 0. 

The main goal of this Appendix is to identify the terms in (0), of order e 4 and higher, 
that make the dynamics generic. In particular we find those terms that contribute to the 
stability of the splay state. We shall then make some comments concerning the incoherent 
manifold. 

B.l Floquet exponents of the splay state 

The splay state is characterized by 

_ / 1 ]£m= ...-2N,-N,0,N,2N,... , , 

Zm ~ \ otherwise 

From the product rule, we see that terms in Eq. fl55|) with more than one factor of z m or 
with m 7^ pN do not contribute to its stability. The terms with a single z m factor have the 
form e" J z m . The lowest order term with no factors of z m with m ^ pN is s 2N e- im iz N . 
(The \zn\ 2 term has order e 2N+2 , since there must be a factor of |ee J '| 2 .) Hence, the stability 
of the splay state is determined by the terms 

6j = c + (s 2 c l e~^z l + e%e- 2W >z 2 + ■■■ + e 2N c N e- N ^z N + c.c.) + 0{e 2N+2 ) (58) 

The terms listed can all be written in as a "pairwise coupled" system: 

1 N 

e i = JfY, G V*- e i)> ( 59 ) 

k=l 

where the real-valued, 27r-periodic coupling function is 

N N 

G{6) = Yl £ 2n c n e ine = c + ]T e 2n {a n cos(n9) + b n sm(n9)} (60) 

n=— JV n=l 



where c_ n = c n , a n = 2Re(c n ), and b n = — 2Im(c n ). See equation (p3|) , where the terms of 
order e A in the pairwise coupled system are given. 
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We emphasize that the averaged equations are not pairwise coupled in the Josephson 
junction system. For example, at order e A there are terms like e~ 2%di z\ and e~ ldi Zi\zi\ 2 on 
the right hand side of ( pH ) that cannot be written as pairwise coupling. But if we truncate 
terms of order e 2N+2 , then the stability of the splay state is determined by pairwise coupled 
terms. In fact, we will find that the the splay state is typically hyperbolic if we truncate 
terms of order e N+1 . 

Note that the Fourier coefficients c n in ([58]) are even functions of e. For example, C\ = 
a + e 2 /3 + 0(e 4 ), where a and (3 are defined in Eq. (Q). In the perturbation expansion, we 
did not compute the (3 correction, and we use c\ to stand for the constant a in (Q). 

In this Appendix we concentrate on the Floquet exponents rather than the multipliers. 
Due to the phase shift symmetry of the averaged equations, the Floquet exponents are simply 
the eigenvalues of a constant matrix: hence we call them A. The Floquet multipliers are 
fi = exp(AT), where T = 2-k/cq is the period of the oscillation. 

Ashwin and Swift ||AshfcSwi92|| have computed the Floquet exponents of the splay solu- 
tions to pairwise systems. Fortunately, this result can be used to compute the approximate 
stability of the splay when system (^) is truncated to order e 2N . The Floquet exponents of 



the splay state in (|59D are 

1 N 

A m = -E G'(2nk/N)(e l2 * km / N - 1), (61) 

k=l 

where m is any integer satisfying —N/2 < m < N/2. There are exactly N Floquet exponents; 
Ao = always, and if N is even then A7V/2 is real. The rest of the exponents come in 
complex conjugate pairs: A_ m = A m for < m < N/2. Thus we only need consider 
1 < m < floor [N/2]. The exponent A m corresponds to the perturbation of z m away from 
zero. 

It is convenient to cast these results in terms of the Fourier coefficients of the coupling 
function G, since these are the numbers that we get from the perturbation expansion. The 
Floquet exponents of the splay solution are 

1 N 00 

A m = -J2 J2 inc n e 2n e i2nnk / N (e l2 * mk / N - 1) (62) 

N k=l n=-oo 

= -t(j2 n£2nC -n + T, n£2nC n) (63) 
\n=m n=0 / 

The sum over n = m means that all n congruent to m modulo iV contribute to the sum. 
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To compute the first nonzero contribution to the Floquet exponents the epsilon expansion 
should be carried out to order iV (if N is even) or N — 1 (if N is odd). Then truncate so 
that c n = if n > N/2, and the Floquet exponents of the splay solution are 

A m = —ime 2m c_ m if \m\ < N/2, and , , 

\ N/2 = (N/2)e N b N/2 if N is even. ^ ' 

The stability of the splay is determined by the real part of the exponent, or equivalently the 

magnitude of the Floquet multiplier: 

|/i± m | = exp(s 2m m7ib m /co) if \m\ < N/2, and , , 

= e *-p( eN NTtb N / 2 /c ) if N is even. 

For N = A oscillators, the critical multiplier fi 2 is determined by the single coefficient b 2 = 

— 2Re(c2), so most of the terms obtained in the averaged equation ( [55] ) at order e A can be 

ignored, as we did in Section 4 by setting Ri = (which is the same as setting z\ — 0). If one 

were brave enough to try to compute c% (for N > 6 oscillators), one should set z\ = z 2 = 

for the calculation, and so on to higher order. 

B.2 The incoherent manifold and antiphase pairs 

Is the incoherent manifold invariant? We show that the answer is no for N > 5. But the set 
of "antiphase pairs," that coincides with the incoherent manifold when N = 4, is shown to 
be dynamically invariant in the averaged equations (|55|) . 

Recall that the incoherent manifold is the set Z\ — 0. The manifold is dynamically 
invariant if ii| 21 =o = 0. The phase shift symmetry and the e ordering of equation (|55|) imply 
that k\ is a linear combination of invariant functions times the terms: 

2 — 4 — 2 4 — 3 4 — frr\ 

z 1 , e z 1 z 2 , e z 1 z 3 , e z 2 z x , e z 2 z 3 , ■■■ (66) 

We see that the incoherent manifold is not dynamically invariant in general due to terms like 
Z2Z3, which contain no factors of z\. The term can be traced to the term proportional 
to v in equation (|5"3]). 

There is, however, a dynamically invariant set. When iV is even, define 

V = • • • , On) : 9j can be permuted so that 9j + N/2 = @j + 71 f° r J = 1,2,..., N/2} 

(67) 

In other words, V has N /2 antiphase pairs. It is easy to see that z m = if m is odd for any 
state in V , because 



37 



System fl55|) can be transformed to an infinite system of ODEs for the moments z m . The 
phase shift invar iance implies that any polynomial term on the right hand side of z m , where 
m is odd, has at least one factor of z p , where p is odd. Hence z m \-p = for odd m, and the 



set V is dynamically invariant in the averaged equation (55). We expect a corresponding 
invariant set in the full Josephson junction equations, but this is not a rigorous proof, since 
it relies on the phase shift symmetry. (See sec. 2.1 of ||AshfeSwi92|| .) 



For N = 2 or iV = 4, the incoherent manifold is V, so it is invariant in ([55]). The 
incoherent manifold is also invariant when N = 3, since it contains only the splay solution. 
For N > 5, the invariant manifold z\ = is not dynamically invariant in (^), but since it 
is invariant and normally hyperbolic in the truncation to "the" averaged equation, there is 
a nearby invariant set when e is sufficiently small. 

B.3 The averaged equations when (3 = 

As we have discussed, the Josephson junction equations with McCumber number j3 = 
have a rather degenerate structure. In particular, the incoherent manifold is invariant. We 
conjecture that the averaged equations in this case are "centroid coupled:" 

oo 

Q 3 = ]T e 2n c n e- in ^zl (69) 

n=— oo 

where c_ n = are different coefficients from those in the pairwise coupled system (|59]). We 
call this a centroid coupled system because z\ is the centroid of the N oscillators placed on 
the unit circle. None of the higher moments, ^2,^3, etc., appear. For a centroid coupled 
system the incoherent manifold, defined by z\ = 0, is dynamically invariant and foliated by 
periodic solutions, all with the same period. There is no "drift" on the incoherent manifold; 
in other words all the moments z m are constant when Z\ = 0, and the Floquet exponents of 
the splay state are all 0, except for a single pair giving the stability normal to the incoherent 
manifold. 

A centroid coupled system can have bistability (or bi-instability) between the in-phase 
solution and the incoherent manifold. For example the e~ 2iej zf term contributes to the 
Floquet exponents of the in-phase solution but not the splay solution. Such bistability is 
observed in a Josephson junction series array with (3 = [TsafcSch92|| . Note that a pairwise 
coupled system cannot describe the array in the limit, since a pairwise coupled system with 
a foliated invariant manifold must be "the" averaged equation (i.e. c n = if |n| > 1 in 
equation (|69|)), which has no bistability. 
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